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Abstract 

Thanks to an algebraic duality property of reduced states, the 
Schmidt best approximation theorems have important corollaries in 
the rigorous theory of two-electron moleculae. In turn, the "har- 
monium model" or "Moshinsky atom" constitutes a non-trivial lab- 
oratory bench for energy functionals proposed over the years (1964- 
today) , purporting to recover the full ground state of the system from 
knowledge of the reduced 1-body matrix. That model is usually re- 
garded as solvable; however, some important aspects of it, in particu- 
lar the exact energy and full state functionals — unraveling the "phase 
dilemma" for the system — had not been calculated heretofore. The 
solution is given here, made plain by working with Wigner quasiprob- 
abilities on phase space. It allows in principle for thorough discussions 
of the merits of several approximate functionals popular in the theo- 
retical chemical physics literature; in this respect, at the end we focus 
on Gill's "Wigner intracule" method for the correlation energy. 



1 Introduction 

At the heart of theoretical chemistry sits the Coulson program of replac- 
ing the wave function by the (reduced) 12-variable 2-body matrix p2 [1]. 
For this function of twelve variables the basic energy functional is trivially 
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known. Unfortunately, the A^-representability problem for p2 has never been 
solved in an usable way. Density functional theory (DFT) was unleashed 
by the Hohenberg-Kohn discovery that somehow the same information as in 
the wave function or in the 2-matrix is stored in the humble 3-variable elec- 
tronic density. However, the energy functional encoding such information is 
unknown, and must have strangely non-local properties — in relation with 
the existence of negative ions, for instance [2]. A natural alternative is the 
Gilbert program [3] based on the 6- variable 1-body matrix pi, for which the 
A^-represent ability question is solved.^ 

Now, possession of the 1-body matrix for an A"-electronic system does 
not allow exact inference of the corresponding 2-body matrix p2, or of the 
correlation energy. It is natural to diagonalize pi and seek to expand p2 in 
terms of eigenfunctions of pi ("natural orbitals") and its eigenvalues < 
Vj < 1 ("occupation numbers"), with 'Ylij^j = ^- Over the years, starting 
with the work by Miiller [8], approximate relations between p2 and p\ and 
energy functionals based on this spectral analysis of pi have been suggested 
and tried, with various results. 

Two-electron systems constitute an exception to our ignorance. For them 
it is feasible to a large extent to reconstruct p2 from pi, and so fulfil the 
Coulson and Gilbert programs. As we shall recall, that it may be so stems 
simply from Schmidt's approximation theorems [9]. For a reduced 1-density 
of the kind 

Pi{x,x') = {^i^v+iiii)pi{,r,r') = {^i\v + iiiv)^n,(Pj{r)(P]{r'), 

^After a long period of stasis, recently there has been progress on the represent ability 
problem for 2-matrices, and on pending questions of represent ability for 1-matrices [4-7] . 
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with rij — 1, the corresponding 2-density matrix is given by 



P2{Xi, X2, X[, X'2) = 2 (^1^2 - 1^2) (tl42' - il'W 




X (Mfi) Mfi) Mfi) 



) 



(1) 



V 



X {<p 



) 



V 



Alas, the recipe, although exact, is underdetermined: of the Cj we only know 
that |cjp = rij. This is a "phase dilemma" of 1-body matrix functional 
theory. We work here only with states described by real wave functions: this 
still leaves us with an infinite number of signs to account for. 

The harmonium model originally proposed by Moshinsky [10] consists of 
two spin-i fermions trapped in a harmonic potential and repelling each other 
with a Hooke's law force, as well. It constitutes an analogue to a two-electron 
atom, helpful to illustrate the main ideas of reduced density matrices and 
correlation energy in an exactly solvable context. Its ground state, as well as 
the reduced density matrix and the pair distribution, exhibiting correlation, 
are well known in the standard wave function formalism [11]. Recently, this 
model has been recruited for suggesting approximate Ansatze for correlation 
energy density [12], and other purposes [13]. 

Nevertheless, within the standard formalism it is not at all obvious how 
to go about the sign conundrum. We manage here to sidestep all the difficul- 
ties, and exemplify the two-electron state theorem for harmonium, by work- 
ing with Wigner functions on phase space instead. Mathematically, Wigner 
phase space states are equivalent to density matrices, sporting of coiuse the 
same number of variables. However, our approach calls for a different type of 
physical intuition, encapsulating some of the non-locality which is the bane 
of DFT. 

Like one-body density matrix functional theory, Wigner functional the- 
ory sits midway between ordinary DFT and the Coulson proposal. Currently 
Gill and coworkers [14] pursue a vigorous program to understand correlation 
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energies via Wigner intracules: our method also sits transversally midway 
between the 6-variable Gilbert and Gill programs. We mention as well a 
different fusion of the Schmidt and Wigner ways, related to information the- 
ory [15], and a recent quantum phase space view of harmonium by Dahl [16]. 

The plan of the article: we start by recalling the classical contribution by 
Schmidt, with a pedagogical bent. We show how to derive from it the exact 
energy functional for two-electron systems, determined up to an unknown 
configuration of signs. After that, we describe harmonium by phase space 
methods; these provide the fast lane towards an elegant and complete de- 
scription of this system. Finally, we examine Gill's Wigner intracule method 
for correlation energy from our results on harmonium. 

2 Schmidt's theorem and its meaning for hehumahkes 

Little more than singular value decomposition of matrices is involved here. 
Let 'i/{xi,X2) be the wave function describing a pure state of any system 
with two electrons — including the helioid series, starting with H~,He, and 
molecules like H2. For simplicity, and because in absence of magnetic fields 
this is actually the case for such systems, we assume that ^ is real. Let us 
expand ^ in suitable orthonormal basis (ONB), also taken real, 

^{xi,X2) = I '^)Mxi)hj{x2) 

ij 

■= ih{xi) f2{xi) fsixi) 



Suppose the eigenvalue problem for the (real positive symmetric) matrix CC* 
is solved by 

CC^Uj = ajUj. 

We order the singular values: cii > (T2 > era > • • • > 0. Then the Uj are the 
columns of an orthogonal matrix U. Define now Vt = C^Ui/ai. These are at 
least partially the columns of another orthogonal matrix, say V , since 

_ {uj I CC\) _ a%Ui I u^) _ ^ 



=■ 2lfi^^^)'^iohj{x2) 

ij 

^hi{x2)\ 

. „ ^2(3^2) 
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Add to them any orthonormal basis for ker C. Notice moreover that 



u\Cvj = —u\CC%j = SijGi, or U^CV = 



(72 

(73 



V ■7 

We just have reached the singular value decomposition of (real) matrices: 

c = UTy\ (2) 

with E diagonal nonnegative and C/, V orthogonal, whose columns are "sin- 
gular vectors" , respectively left-hand and right-hand — see for instance [17, 
Chap. 5]. Notice the non-unicity of V 

Before proceeding, wc remark that Schmidt in [9] showed that the best 
finite-rank approximation (in the sense of any orthogonally invariant norm) 
to C is given simply by 

where S*^^^ is obtained just by cutting out ak+i = <7k+2 = • • • — 0, and 
U^^\ y^'^) are the matrices obtained by conformally keeping in [/, V the first 
k singular vectors. Most important for us, we note as well that CC* (save 
the conventional factor 2) precisely represents the 1-body density associated 
to 

j *(a;, X2)'^{x', X2) dx2 = fi{x) CijCkj fk{x') = af(pi{x)(j)i{x') 

ijk I 

I 

where we have made the spectral expansion of that density: the i^i are the 
chemists' occupation numbers and the (f)i the natural orbitals. We witness 
here a wonderful convergence of mathematical relevance and physical import. 
Schmidt approximation then implies that a state made out of K{K — l)/2 
configurations can be rewritten with K configurations. 



^This decomposition amounts to the canonical "Schmidt scries" for an operator and 
the "Schmidt decomposition" of tensor products, currently taught to students interested 
in entanglement, decoherence and quantum information theories [18]. 
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Going back to (2) and denoting: 



h2 



(P2 



it is clear that C*C is the matrix representation of pi in the ONB { \hj) } 
and CC* the one in the ONB { }, and that 

*(Xi,X2) = ± V (f>i{Xi)ipi{x2) 

i 

-■ y^^±y/i^i/2(j)i{xi) / dx(j)i{x)'^{x,X2). 

i 

The sign indetermination at this stage cannot be removed. The challenge 
is to identify the right configuration of signs for each system from energy 
considerations. 

Since moreover is skewsymmctric. (fi must be a natural orbital itself. 
For instance, for the singlet ground state of the helioids, 

^g,{Xi,X2) = ;^(ti;2 -ilt2)^(n,?'2), 

if we choose 

hi-i(x) = t0i(r), (l)2i(x) = 4,0i(r), 
then ip2i-i(x) ^ l(l)i(r), ip2i(x) ^ -"tMr). 

For the triplet states we have: 

*ti(2;i,X2) = ^1^2"^ di[yi{ri)uji{r2) - uji{ri)-fi{r2)]; 

i 

*to(2;i,X2) = -^{'tii2 + iit2)^di[li{ri)i^i{r2) - uji{ri)-fi{r2)]; 

* i 

^'t-i(a;i,a;2) = iii2^di[%{ri)ui{r2) - (^i{ri)%{r2)]; 



so, if we choose 



then ipli_^ix) = M(r), (pli{x) = -t^iir); 
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and similarly in the other cases. 

The Schmidt theorem has a far-reaching generalization for best approx- 
imations of A^-body wave functions in terms of orbitals of complementary 
reduced density matrices. A good place to see that in action is [11]. 



3 Harmonium on phase space 

The Hamiltonian for harmonium in Hartree units is 

^= 2"+ 2"+ 2(^1 + ^2) - ^?^i2- (3) 
Introduce extracule and intracule coordinates, respectively given by 

H = (n + r2)/^/2, r = (n - r2)/v^, 
with conjugate momenta 

P = (Pi + P2)/V2, p = (pi - P2)/ V2. 



Therefore 



^ ^ „ P2 ^2^2 2 2^2 



2 2 2 2 

We have introduced the frequencies u = y/k and n = \/k — 5; for both 
particles to remain in the potential well, assume 6 < k. 

The translation to phase space of the analysis in the previous section was 
accomphshed in [19]. One could depart somewhat from the treatment in this 
reference, though. There, following in particular Springborg and Dahl [20], 
by outright transfer from the density matrix formalism, we regard the 1- 
body and 2-body Wigner functions defined below as respectively 2x2,4x4 
matrices. It is however more natural, as advised by Wigner himself [21] and 
effected in [22], to describe spin by combining row and column components 
to have simpler transformation properties; these are none other than the 
ones corresponding to the tensor representations of the rotation group. In 
particular for the reduced 1-density for a two-electron system in a singlet 
spin state, the spin factor: 

^(titi'+;i;i') 
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is just the identity scalar, and can be omitted. In general, for the 1-body 
Wigner function one organizes the spin factors in terms of the Cartesian 
components of the spin operator vector. 

We may concentrate on the spinless part of the ground state henceforth. 
In general, the spinless Wigner quasiprobability (normalized to one) corre- 
sponding to a (real) 2-particle wave function ^ is given by 

P^s.{ri,r2;pi,P2) (4) 
= ^ y" P2(ri - zi, r2 - Z2; ri + Zi, r2 + Z2) e^^^^-^^+P-^^) dh^ d'z2, 

with P2{Xi,r2, 'r'i,r'2) = ^(ri, r2)^(ri, r2). The definition extends to transi- 
tion matrices |$)($'| as well: 

P<E,$/(ri,r2;pi,p2) 
= ^J $(ri - zi, r2 - Z2)$'(ri + Zi, r2 + Z2) c'^'^p'-'^'+p'-'''^ dhi d^Z2. 

For harmonium the Wigner distribution factorizes into extracule and intra- 
cule parts: 

Pgs{ri,r2;PuP2) = ^^xpl — — 1 expl — — 1. (5) 

This is reached most efficiently and elegantly by the methods of phase space 
quantum mechanics [23] . One can now obtain p2 by the inverse formula of (4) 
— easily provided for by Fourier analysis. The pairs density P2(?"i, ^2, fi, ^2) 
is recovered by integration over the momenta. 

The reduced 1-body phase space (spinless) quasidensity for the ground 
state dgs is obtained, as in the standard formalism, by integrating out one 
set of variables. 



^gs(r-;p) = ^(^ 



(^ + /^)V 



It is a simple exercise to recover from this Eq. (2-68) for pi{ri,r'i) in the 
treatment in [11]. The marginals of dgg give the electronic density and mo- 
mentum density. 

It should be recognized that, while Pgs is a pure state, mathematically 
dgs describes a mixed state. For Gaussians on phase space, such as Pgs 
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and dgs too, there are simple rules to determine whether they represent a 
pure state [24], a mixed state [25], or neither. Writing q = (r"i,r2), tt = 
(pi,P2), u = {q,n), we find Pgs(u) = Tr-^e""-^" = Tr-^e"'?-^'?-"-^"''^ where, 
amusingly, 

^-^v A-^=u^";+^"; "^"'rO. (7) 

We see that the matrix F corresponding to formula (5) is symmetric and 
symplectic, and therefore represents indeed a pure state. This is not the case 
for dgs- Thus our purpose in recovering Pgg from knowledge of dgs alone is 
akin to putting Humpty Dumpty together again! 

4 ... who traced these signs? 

Since all the relevant quantities factorize, in this section we work in one 
dimension for notational simplicity. Write u = {r,p). For the ground state 
on phase space formula (1) goes [19] into: 

n / .V (tl4'2 ~ 4'lt2) (tl42' ~ 4'l't2') 

-PgsK,ii2;spm) = 



X ^ CrCsfrs{ui)frs{u2), (8) 
r,s=0 

where CrCg = {±^/n^){±^/n^), with the signs to be determined. The frs 
are Wigner eigentransitions, the functions on phase space corresponding to 
transitions between states. A similar series for triplet states is established 
and investigated in [26]. 

The explicit form of the frs will turn out to be well known. But first we 
have to find the good coordinates. According to a theorem in [25], the real 
quadratic form in the exponent of dgs must be symplectically congruent to a 
diagonal one. Wc perform the transformation 

{Q,P) '■— or, in shorthand, U — Su, 

where S is symplectic. Introducing as well the parameter A := 2y/ujfl/ (uu+fj,), 
the 1-quasidensity takes the simple form 

dgsiU) := dgs{u{U)) = —e-^''\ 

TT 
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Prom now on, we assume in the notation S > 0; the treatment for the "at- 
tractive" case u < fi foUows parallel lines. Under that proviso, we may also 
write A =: tanh(/3/2), so that 

P = log r- = 2 log — and smh /3/2 = = . 

Prom the series formula. 



n=0 

it follows that 



d,,{U) = - sinh|5^(-l)'-L,(2t/^)e-^'e-(^'^+^)'^/l 

One recognizes the basis of Wigner eigenfunctions on phase space correspond- 
ing to oscillator states [23] : 

frr{U) = -{-iyLr{2U')e-''"; (9) 

TT 

thus we realize that dgs is in thin disguise a Gibbs state [25] with inverse 
temperature (3: 

d^{U) = 4sinh I Yl e-^'^-^'^^/'frriU). (10) 

r=0 

With this we have identified the natural orbitals in the U variables. Their 
occupation numbers are 

_ 2r+l 

n, = 2sinh^e-(2'-+^)^/2 



Ay/uJJI /ycJ- 

U - H \^/uJ + ^/JiJ 



Notice that no = 1 — e^^ = Z^^{I3), for Z the system's partition function; 
also ^^rir = (1 — e~^) e~^^ = 1. These have nice square roots: 



10 



Finally, for r > s one has [23] : 

frs(U) := l(-l)^/4 (2U^f-'^/^e-'^'-'^^Ll-'{2U^)e-^\ 



TT V r! 



where := aictan{P/Q). The fsr are complex conjugates of frs- 

We turn now to the energy functionals, which should allow us to determine 
the correct signs in (8). Of course, still working in dimension one, the total 
energy in our units is just one half of the sum of the frequencies; and this 
is confirmed by the well-known formula for the Gibbs ensemble (consult [27] 
for the phase space derivation) represented by (10). To wit, 

Eg, = E[d^,] = 2^(^^ + + ^ 



Ve/^-1 2 

Prom the viewpoint of DFT, the most interesting part of the energy is 
the electron interaction £^2gs- The 1-body Hamiltonian is given by 

It is a simple exercise to obtain the 1-body energy i^igs by integrating ex- 
pression (6) with this observable: 

^'^^ - 2 + 

It is easily checked that this equals 2 nrE„-, where denotes the 1-body 
energy associated to each natural orbital in (9). Thus it must be that 

F -F F - 

-C'2gs — -C/gs - -C/lgs — ■ 

Now we go for the two-body contributions of the natural orbitals. The 
interelectronic repulsion potential in (3) is (/x^ — a;^)r^2/4, so they are of the 
form Ylirs^TCsLsri with the Lgr given by: 

,2 , ,2 



Lsr = — ^ — / fsr{qi;pi)fsr{q2;p2){qi - q2f dqi dq-2 dpi dp2 

hs{Ql)hr{Qi){Qi - Q2fK{Q2)hr{Q2) dQi dQ2. 
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Here hj. are the usual orthogonal harmonic oscillator eigenfunctions for unit 
frequency. We consider the diagonal r — s first, whereby 



1^ 



and thus 



16a;// 



Apt 



2uj 



We have used that the expected value of for a harmonic oscillator eigen- 

state isr+l and that the expected value of Q is zero. Notice that — — — < 1. 

2a; 

We have been considering /j, < to. For attraction between the electrons, this 
contribution changes sign. 

A non- vanishing contribution of the off-diagonal part may then come only 
from the terms 



CO 



/x^ 



hr{Q)hr+i{Q)Q dQ 



We compute: 



a;2-//2 



a;2 



OO r » 

y/nrUr+l j hr{Q)hr+i{Q)Q dQ 

r=0 



a;// {^^^Y 



r=0 



= (a;^ 



U! — fjL U! — fjL 



[^/u+y/|I)^ 16a;/t 



4// 



2a; 



Here we employ ^^o(^ + l)-^*^ — ~ x)~'^. The factor (r + l)/2 comes 
from the definition of the emission operators — {Q — iP)/V2 (or the 
absorption operators), with a^hr — \/r~+T/ij.+i. There is also an overall 
factor of 2 coming from two subdiagonals for each r. Obviously there are no 
other contributions. 

In conclusion, to minimize the energy we now have to choose minus signs 
whenever s — r±l. This conclusion does not depend on whether the electrons 
attract or repel each other. It determines all the signs in (8). Then the total 
2-body energy comes out as 



/t^ — a;^ a; + /t 



a;^ — /t^ a; — /i 



12^ -U^ 



4// 



2a; 



4// 



2a; 



4// 
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as it ought to be. 

Theorem 1. The ground state of harmonium is equal to 

Pgs{Ul,U2) = (-)''^V"'r^s frs{Ui)frs{u2), 

where the frs are the oscillator eigentransitions up to a symplectic change of 
basis, and the Ur have been given in (11). 

According to the spirit of DFT, the above result has been derived from 
knowledge of dgg alone; the advantage of working with Wigner distributions 
being that the necessary transformations arc completely transparent. How- 
ever, since we had Pgs explicitly beforehand, it is nice to make the check. 
In (8) sum over each subdiagonal, where r — s = / > 0: 



1 ^ 

1 — s=l 



^Z. (/ + ,)! ^ ^sWJ^.l^^2J- ^2 e '^^sinh(/3/2)y' 
where li is the modified Bessel function. The generating function identity 



h{z)^2Y,h{z) cos(/^) = e^^°^^ 
where the sign rule dictates ^ = ■j?! + ■j?2 + tt, yields for the total sum: 

^-2g-[(!72+[7|)/A+2C/i!72 csch(/3/2) cos(i?i+i92)] (^2) 
^ ^-2g-i[(r2+r2)(a;+M)+{pf+p2)(a;-i+;,-i)]g-rir2(a;-M)gPiP2(M-'-'^-')^ 

which is the correct result. To our knowledge this the first analytical verifi- 
cation of the two-electron state theorem. 

Once the correct sign configuration in the functional is obtained, one 
may compare it with standard approximate workhorses like the Miiller func- 
tional [8,28]. Miiller's is a clever variation on the standard Hartree-Fock 
functional. For a two-electron system in the phase space approach, it adopts 
instead of (8) the form: 

Pm{ux, U2\ spin) = ^ (tiia - ^vii' - iv^i^ 

(oo oo \ 

^ nrnsfrr{Ul)fss{u2) - \/n^frs{Ul)fsr{u2)y 
r,s=0 r,s=0 ^ 
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It is well known that for Coulombian systems the Miiller functional tends to 
overbind. The rigorous proof of this fact for real two-electron atoms given 
in [19] is much more transparent than the one in [28]; it shows that definite 
positivity of the Coulomb potential does play a decisive role, whereas the 
extra minus signs in Miiller's functional do not. For these very reasons the 
Miiller functional's tendency to overcorrelate needs reexamination in harmo- 
nium. This is taken up in [26]. 

5 On the correlation functional 

The Wigner quasiprobability for the best Hartree-Fock state for harmonium 
is given by 

^ 1. g-(i?2+r2)^(a;2+Ai2)/2g-(p2+p2)/^(a,2+^2)/2 

We refer to [10] for the derivation of the Hartree-Fock wave function for 
harmonium, from which the above follows by (4). From the above, comparing 
with the exact state under the form (5) or (12), a "Moshinsky hole" in phase 
space is clearly visible. 

The correlation energy Ec is defined as the difference between the true en- 
ergy of the electronic system and the Hartree-Fock energy. For harmonium, 
use of Phf gives 



2 2 2 

^ g-(H2+,2)^(^,+^.)/2 ^-iP^+p-y,/i^:pTJ^ d'Pd'Rd'pd'r 



4 4 + /x2)/2 4 4,/{u^ + ^^)/2 

and so the correlation energy is 



■— ^0 ~ ^HF — <J 



32a;3' 
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as (5 4, keeping oo fixed, or as w f oo keeping 5 fixed. If we understand 
the latter as the "high-density" fimit, the correlation energy vanishes in it. 
In [29] the high energy limit is understood as the limit /i t c>o, with u kept 
constant. In both cases the behaviour of correlation is at variance with that 
of the Coulomb potential [30] . 

As indicated in the introduction, Gill's program for the correlation en- 
ergy [14,31] is based exclusively on the Wigner intracule. Since relative mo- 
mentum is as important as relative position in determining interelectronic 
correlation, it is certainly appealing to study correlation on phase space. 
More precisely, the full Wigner distribution is reduced to the relative vari- 
ables ri2 = Ti — r2, Pi2 —Pi—P2- Then it is conjectured that the correlation 
energy is a universal functional of the "Omega intracule" fi, that is, the in- 
tracule on phase space after reduction to the variables ri2,pi2, w, with w the 
angle between ri2 and pi2, of the form: 

Jo Jo Jo 
or even a functional of the lower intracules: 

poc poo poo 

Ec^ / ^iri2,Pi2)G{ri2,Pi2)dpi2dru; Ec= A{s)G{s)ds, 
Jo Jo Jo 

with s :— ri2Pi2- In our case, recalling r — Tyil^J^-, P — Pi2/V^, the exact 
intracule takes the form 

i_e-M.^-pVM _^ 1^2 2 ^-^rl,l2^-,l,l2^ ^. W,^{n2,Pl2), 

after integration over the angles. While the Hartree-Fock intracule is: 

TT 

which, as observed in [29] , is of the same form as the exact one when replac- 
ing 5 by S/2. As a consequence the "action intracules": 

A{s) ^ - r rl, 4- e-'--/^e-V2..?. ±11 = —K,{s) 

TTJo ^2 ^2 TT 

for any value of •, are precisely the same for the exact and Hartree-Fock 
cases. Note that in the exact solution only the intracule part of the Wigner 
quasiprobability changes in harmonium, while the Hartree-Fock method mod- 
ifies both the extracule and the intracule; and then both contribute to the 
correlation energy. However, in Gill's I^-functions the footprints of both 
harmonium frequencies are jointly kept. 
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6 Conclusion and outlook 



We have beaten a path from Schmidt series to the exact energy functional 
for two-electron systems, and fully determined the latter for the harmonium 
model of an atom, by working with Wigner distributions on phase space. We 
hope this paper helps to understand why this kind of approach is bound to 
grow in importance in the coming years. 
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